Plots
rm(list=ls());gc()## used (Mb) gc trigger (Mb) max used (Mb)
## Ncells 426256 22.8 893647 47.8 648454 34.7
## Vcells 764472 5.9 8388608 64.0 1661560 12.7load("G:/Mi unidad/paperestallido/Definitive_models_2021.RData")
#xaringan::inf_mr()
#setwd("C:/Users/CISS Fondecyt/Pictures")
#C:\Users\CISS Fondecyt\Pictures#arriba puse algunas opciones para que por defecto escondiera el código
#también cargue algunos estilo .css para que el texto me apareciera justificado, entre otras cosas.
local({r <- getOption("repos")
r["CRAN"] <- "http://cran.r-project.org"
options(repos=r)
})
`%>%` <- magrittr::`%>%`
copy_names <- function(x,row.names=FALSE,col.names=TRUE,dec=",",...) {
if(class(ungroup(x))[1]=="tbl_df"){
if(options()$OutDec=="."){
options(OutDec = dec)
write.table(format(data.frame(x)),"clipboard",sep="\t",row.names=FALSE,col.names=col.names,...)
options(OutDec = ".")
return(x)
} else {
options(OutDec = ",")
write.table(format(data.frame(x)),"clipboard",sep="\t",row.names=FALSE,col.names=col.names,...)
options(OutDec = ",")
return(x)
}
} else {
if(options()$OutDec=="."){
options(OutDec = dec)
write.table(format(x),"clipboard",sep="\t",row.names=FALSE,col.names=col.names,...)
options(OutDec = ".")
return(x)
} else {
options(OutDec = ",")
write.table(format(x),"clipboard",sep="\t",row.names=FALSE,col.names=col.names,...)
options(OutDec = ",")
return(x)
}
}
}
unlink('*_cache', recursive = TRUE)
if(!require(pacman)){install.packages("pacman")}
pacman::p_unlock(lib.loc = .libPaths()) #para no tener problemas reinstalando paquetes
knitr::opts_chunk$set(
echo = TRUE,
message = FALSE,
warning = FALSE
)
#dejo los paquetes estadÃsticos que voy a utilizar
if(!require(plotly)){install.packages("plotly")}
if(!require(lubridate)){install.packages("lubridate")}
if(!require(htmlwidgets)){install.packages("htmlwidgets")}
if(!require(tidyverse)){install.packages("tidyverse")}
if(!require(gganimate)){install.packages("gganimate")}
if(!require(readr)){install.packages("readr")}
if(!require(stringr)){install.packages("stringr")}
if(!require(data.table)){install.packages("data.table")}
if(!require(DT)){install.packages("DT")}
if(!require(ggplot2)){install.packages("ggplot2")}
if(!require(lattice)){install.packages("lattice")}
if(!require(forecast)){install.packages("forecast")}
if(!require(zoo)){install.packages("zoo")}
if(!require(panelView)){install.packages("panelView")}
if(!require(janitor)){install.packages("janitor")}
if(!require(rjson)){install.packages("rjson")}
if(!require(estimatr)){install.packages("estimatr")}
if(!require(CausalImpact)){install.packages("CausalImpact")}
if(!require(textreg)){install.packages("textreg")}
if(!require(sjPlot)){install.packages("sjPlot")}
if(!require(foreign)){install.packages("foreign")}
if(!require(tsModel)){install.packages("tsModel")}
if(!require(lmtest)){install.packages("lmtest")}
if(!require(Epi)){install.packages("Epi")}
if(!require(splines)){install.packages("splines")}
if(!require(vcd)){install.packages("vcd")}
if(!require(astsa)){install.packages("astsa")}
if(!require(forecast)){install.packages("forecast")}
if(!require(MASS)){install.packages("MASS")}
if(!require(ggsci)){install.packages("ggsci")}
if(!require(Hmisc)){install.packages("Hmisc")}
if(!require(compareGroups)){install.packages("compareGroups")}
if(!require(dplyr)){install.packages("dplyr")}
if(!require(ggforce)){install.packages("ggforce")}
if(!require(imputeTS)){install.packages("imputeTS")}
if(!require(doParallel)){install.packages("doParallel")}
if(!require(SCtools)){install.packages("SCtools")}
if(!require(MSCMT)){install.packages("MSCMT")}
# Calculate the number of cores
no_cores <- detectCores() - 1
cl<-makeCluster(no_cores)
registerDoParallel(cl)
Sys.setlocale(category = "LC_ALL", locale = "english")## [1] "LC_COLLATE=English_United States.1252;LC_CTYPE=English_United States.1252;LC_MONETARY=English_United States.1252;LC_NUMERIC=C;LC_TIME=English_United States.1252"0.1 Final Plots
library(cowplot)
#horizontal
height_y<-16
height_x<-16
size_title<-18
line_size<-1.2 #.8
angle_x<-45
Sys.setlocale(category = "LC_ALL", locale = "english")## [1] "LC_COLLATE=English_United States.1252;LC_CTYPE=English_United States.1252;LC_MONETARY=English_United States.1252;LC_NUMERIC=C;LC_TIME=English_United States.1252"hosp_trauma_plot_red<-
plot(impact3d_hosp_trauma, "original")$data %>%
dplyr::left_join(data15a64_rn[,c("date","rn")], by=c("time"="rn")) %>%
ggplot(aes(x=date))+
geom_line(aes(y=mean), size=line_size, linetype="longdash", color="red") +
geom_line(aes(y=response), size=line_size, color="black")+
geom_ribbon(aes(ymin=lower,ymax=upper),fill="steelblue", alpha = 0.35)+
theme_sjplot()+
scale_x_date(date_breaks = "3 months",date_labels="%b %y",limits=c(as.Date("2015-01-05"),as.Date("2019-12-23")))+
ylim(c(0,125))+
geom_vline(xintercept = length(data15a64_rn$did[data15a64_rn$did==0])+1, col = "gray50", lty = 2, size=1)+
theme(axis.text.x = element_text(angle = angle_x, vjust = 0.5, hjust=.5,size=height_x),
# theme(axis.text.x = element_blank(),#element_text(angle = 90, vjust = 0.5, hjust=.5,size=9),
axis.text.y = ggtext::element_markdown(size=height_y),
axis.title.y = element_blank(),
axis.title.x = element_blank(),
plot.title = ggtext::element_markdown(size=size_title),
plot.caption = ggtext::element_markdown(hjust = 0, face= "italic",size=9))+
ggtitle("c) Trauma Hospitalizations")
hosp_resp_plot_red<-
plot(impact3d1_hosp_resp, "original")$data %>%
dplyr::left_join(data15a64_rn[,c("date","rn")], by=c("time"="rn")) %>%
ggplot(aes(x=date))+
geom_line(aes(y=mean), size=line_size, linetype="longdash", color="red") +
geom_line(aes(y=response), size=line_size, color="black")+
geom_ribbon(aes(ymin=lower,ymax=upper),fill="steelblue", alpha = 0.35)+
theme_sjplot()+
scale_x_date(date_breaks = "3 months",date_labels="%b %y",limits=c(as.Date("2015-01-05"),as.Date("2019-12-23")))+
ylim(c(0,125))+
geom_vline(xintercept = length(data15a64_rn$did[data15a64_rn$did==0])+1, col = "gray50", lty = 2, size=1)+
theme(axis.text.x = element_text(angle = angle_x, vjust = 0.5, hjust=.5,size=height_x),
axis.text.y = ggtext::element_markdown(size=height_y),
axis.title.y = element_blank(),
axis.title.x = element_blank(),
plot.title = ggtext::element_markdown(size=size_title),
plot.caption = ggtext::element_markdown(hjust = 0, face= "italic",size=9))+
ggtitle("d) Respiratory Hospitalizations")
cons_trauma_plot_red<-
plot(impact3d1_cons_trauma, "original")$data %>%
dplyr::left_join(data15a64_rn[,c("date","rn")], by=c("time"="rn")) %>%
ggplot(aes(x=date))+
geom_line(aes(y=mean), size=line_size, linetype="longdash", color="red") +
geom_line(aes(y=response), size=line_size, color="black")+
geom_ribbon(aes(ymin=lower,ymax=upper),fill="steelblue", alpha = 0.35)+
theme_sjplot()+
scale_x_date(date_breaks = "3 months",date_labels="%b %y",limits=c(as.Date("2015-01-05"),as.Date("2019-12-23")))+
ylim(c(0,1600))+
geom_vline(xintercept = length(data15a64_rn$did[data15a64_rn$did==0])+1, col = "gray50", lty = 2, size=1)+
theme(axis.text.x = element_text(angle = angle_x, vjust = 0.5, hjust=.5,size=height_x),
axis.text.y = element_text(size=height_y),
axis.title.y = element_blank(),
axis.title.x = element_blank(),
plot.title = element_text(size=size_title),
plot.caption = element_text(hjust = 0, face= "italic",size=9))+
ggtitle("a) Trauma Consultations")
cons_resp_plot_red<-
plot(impact3d_cons_resp, "original")$data %>%
dplyr::left_join(data15a64_rn[,c("date","rn")], by=c("time"="rn")) %>%
ggplot(aes(x=date))+
geom_line(aes(y=mean), size=line_size, linetype="longdash", color="red") +
geom_line(aes(y=response), size=line_size, color="black")+
geom_ribbon(aes(ymin=lower,ymax=upper),fill="steelblue", alpha = 0.35)+
theme_sjplot()+
scale_x_date(date_breaks = "3 months",date_labels="%b %y",limits=c(as.Date("2015-01-05"),as.Date("2019-12-23")))+
ylim(c(0,1600))+
geom_vline(xintercept = length(data15a64_rn$did[data15a64_rn$did==0])+1, col = "gray50", lty = 2, size=1)+
theme(axis.text.x = element_text(angle = angle_x, vjust = 0.5, hjust=.5,size=height_x),
axis.text.y = element_text(size=height_y),
axis.title.y = element_blank(),
axis.title.x = element_blank(),
plot.title = element_text(size=size_title),
plot.caption = element_text(hjust = 0, face= "italic",size=9))+
ggtitle("b) Respiratory Consultations")
rate_trauma_plot_red<-
plot(impact3d_ratio, "original")$data %>%
dplyr::left_join(data15a64_rn[,c("date","rn")], by=c("time"="rn")) %>%
ggplot(aes(x=date))+
geom_line(aes(y=mean), size=line_size, linetype="longdash", color="red") +
geom_line(aes(y=response), size=line_size, color="black")+
geom_ribbon(aes(ymin=lower,ymax=upper),fill="steelblue", alpha = 0.35)+
theme_sjplot()+
scale_x_date(date_breaks = "3 months",date_labels="%b %y",limits=c(as.Date("2015-01-05"),as.Date("2019-12-23")))+
ylim(c(0,400))+
geom_vline(xintercept = length(data15a64_rn$did[data15a64_rn$did==0])+1, col = "gray50", lty = 2, size=1)+
theme(axis.text.x = element_text(angle = angle_x, vjust = 0.5, hjust=.5,size=height_x),
axis.text.y = element_text(size=height_y),
axis.title.y = element_blank(),
axis.title.x = element_blank(),
plot.title = element_text(size=size_title),
plot.caption = element_text(hjust = 0, face= "italic",size=9))+
ggtitle("e) Trauma Hospitalizations per 1,000 consultations")
rate_resp_plot_red<-
plot(impact3d_ratio_resp, "original")$data %>%
dplyr::left_join(data15a64_rn[,c("date","rn")], by=c("time"="rn")) %>%
ggplot(aes(x=date))+
geom_line(aes(y=mean), size=line_size, linetype="longdash", color="red") +
geom_line(aes(y=response), size=line_size, color="black")+
geom_ribbon(aes(ymin=lower,ymax=upper),fill="steelblue", alpha = 0.35)+
theme_sjplot()+
scale_x_date(date_breaks = "3 months",date_labels="%b %y",limits=c(as.Date("2015-01-05"),as.Date("2019-12-23")))+
ylim(c(0,400))+
geom_vline(xintercept = length(data15a64_rn$did[data15a64_rn$did==0])+1, col = "gray50", lty = 2, size=1)+
theme(axis.text.x = element_text(angle = angle_x, vjust = 0.5, hjust=.5,size=height_x),
# theme(axis.text.x = element_blank(),#element_text(angle = 90, vjust = 0.5, hjust=.5,size=9),
axis.text.y = element_text(size=height_y),
axis.title.y = element_blank(),
axis.title.x = element_blank(),
plot.title = element_text(size=size_title),
plot.caption = element_text(hjust = 0, face= "italic",size=9))+
ggtitle("f) Respiratory Hospitalizations per 1,000 consultations")
#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:
#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:
plot1b_red<-
plot_grid(
cons_trauma_plot_red+theme(axis.text.x = element_blank()),
cons_resp_plot_red+theme(axis.text.x = element_blank()),
hosp_trauma_plot_red+theme(axis.text.x = element_blank()),
hosp_resp_plot_red+theme(axis.text.x = element_blank()),
rate_trauma_plot_red,
rate_resp_plot_red,#,
nrow = 3, rel_widths = c(1, 1), rel_heights = c(rep(1,2),1.25)
)
ggdraw(plot1b_red)+
theme(plot.background = element_rect(fill=NA, color = NA))
Figure 1. Estimated Trends of the Outcomes
jpeg("_FigS1_final.jpg", height = 14, width = 23, units = 'in', res = 600)
ggdraw(plot1b_red)+
theme(plot.background = element_rect(fill=NA, color = NA))
dev.off()## png
## 2 #add_sub(, "Note. The first panel shows the data and a counterfactual prediction for the post-treatment period (Blue dashed line);\nBlue area= Prediction intervals.", vpadding=grid::unit(0, "lines"), y = .55, x = 0.03, hjust = 0,size=9)
#plot2 <- cowplot::ggdraw(grid.arrange(p14, p21,p212,p32,p42,p52,p57, ncol = 2, nrow = 4)) +
# same plot.background should be in the theme of p1 and p2 as mentioned above
# theme(plot.background = element_rect(fill=NA, color = NA))
cairo_ps("_FigS1_final.eps",
height = 14, width = 23,
fallback_resolution = 600,
family= "Times")
#horizontal=T,
#paper= "letter",
#pagecentre=F)
ggdraw(plot1b_red)+
theme(plot.background = element_rect(fill=NA, color = NA))
dev.off()## png
## 2 ggsave("_FigS1_final.ps", plot =
ggdraw(plot1b_red)+
theme(plot.background = element_rect(fill=NA, color = NA)),
path = NULL,
scale = 1, width = 23, height = 14, units = "in",
device=cairo_ps,
fallback_resolution = 600,
limitsize = F)#horizontal
height_y<-16
height_x<-16
size_title<-18
line_size<-.8
angle_x<-45
hosp_trauma_plot2_red<-
plot(impact3d_hosp_trauma, "pointwise")$data %>%
dplyr::left_join(data15a64_rn[,c("date","rn")], by=c("time"="rn")) %>%
ggplot(aes(x=date))+
geom_line(aes(y=mean), size=line_size, linetype="longdash", color="darkblue") +
geom_line(aes(y=response), size=line_size, color="black")+
geom_ribbon(aes(ymin=lower,ymax=upper),fill="steelblue", alpha = 0.35)+
theme_sjplot()+
geom_hline(yintercept = 0, col = "black", lty = 3, size=1)+
geom_vline(xintercept = length(data15a64_rn$did[data15a64_rn$did==0])+1, col = "gray50", lty = 2, size=1)+
theme(axis.text.x = element_blank(),
# theme(axis.text.x = element_blank(),#element_text(angle = 90, vjust = 0.5, hjust=.5,size=9),
axis.text.y = element_text(size=height_y),
axis.title.y = element_blank(),
axis.title.x = element_blank(),
plot.title = element_text(size=size_title),
plot.caption = element_text(hjust = 0, face= "italic",size=9))+
ggtitle("c) Trauma Hospitalizations")+
scale_x_date(date_breaks = "3 months",date_labels="%b %y",limits=c(as.Date("2015-01-05"),as.Date("2019-12-23")))
hosp_resp_plot2_red<-
plot(impact3d1_hosp_resp, "pointwise")$data %>%
dplyr::left_join(data15a64_rn[,c("date","rn")], by=c("time"="rn")) %>%
ggplot(aes(x=date))+
geom_line(aes(y=mean), size=line_size, linetype="longdash", color="darkblue") +
geom_line(aes(y=response), size=line_size, color="black")+
geom_ribbon(aes(ymin=lower,ymax=upper),fill="steelblue", alpha = 0.35)+
theme_sjplot()+
geom_hline(yintercept = 0, col = "black", lty = 3, size=1)+
geom_vline(xintercept = length(data15a64_rn$did[data15a64_rn$did==0])+1, col = "gray50", lty = 2, size=1)+
theme(axis.text.x = element_blank(),
# theme(axis.text.x = element_blank(),#element_text(angle = 90, vjust = 0.5, hjust=.5,size=9),
axis.text.y = element_text(size=height_y),
axis.title.y = element_blank(),
axis.title.x = element_blank(),
plot.title = element_text(size=size_title),
plot.caption = element_text(hjust = 0, face= "italic",size=9))+
ggtitle("d) Respiratory Hospitalizations")+
scale_x_date(date_breaks = "3 months",date_labels="%b %y",limits=c(as.Date("2015-01-05"),as.Date("2019-12-23")))
cons_trauma_plot2_red<-
plot(impact3d1_cons_trauma, "pointwise")$data %>%
dplyr::left_join(data15a64_rn[,c("date","rn")], by=c("time"="rn")) %>%
ggplot(aes(x=date))+
geom_line(aes(y=mean), size=line_size, linetype="longdash", color="darkblue") +
geom_line(aes(y=response), size=line_size, color="black")+
geom_ribbon(aes(ymin=lower,ymax=upper),fill="steelblue", alpha = 0.35)+
theme_sjplot()+
geom_hline(yintercept = 0, col = "black", lty = 3, size=1)+
geom_vline(xintercept = length(data15a64_rn$did[data15a64_rn$did==0]), col = "gray50", lty = 2, size=1)+
theme(axis.text.x = element_blank(),
# theme(axis.text.x = element_blank(),#element_text(angle = 90, vjust = 0.5, hjust=.5,size=9),
axis.text.y = element_text(size=height_y),
axis.title.y = element_blank(),
axis.title.x = element_blank(),
plot.title = element_text(size=size_title),
plot.caption = element_text(hjust = 0, face= "italic",size=9))+
ggtitle("a) Trauma Consultations")+
scale_x_date(date_breaks = "3 months",date_labels="%b %y",limits=c(as.Date("2015-01-05"),as.Date("2019-12-23")))
cons_resp_plot2_red<-
plot(impact3d_cons_resp, "pointwise")$data %>%
dplyr::left_join(data15a64_rn[,c("date","rn")], by=c("time"="rn")) %>%
ggplot(aes(x=date))+
geom_line(aes(y=mean), size=line_size, linetype="longdash", color="darkblue") +
geom_line(aes(y=response), size=line_size, color="black")+
geom_ribbon(aes(ymin=lower,ymax=upper),fill="steelblue", alpha = 0.35)+
theme_sjplot()+
geom_hline(yintercept = 0, col = "black", lty = 3, size=1)+
geom_vline(xintercept = length(data15a64_rn$did[data15a64_rn$did==0])+1, col = "gray50", lty = 2, size=1)+
theme(axis.text.x = element_blank(),
# theme(axis.text.x = element_blank(),#element_text(angle = 90, vjust = 0.5, hjust=.5,size=9),
axis.text.y = element_text(size=height_y),
axis.title.y = element_blank(),
axis.title.x = element_blank(),
plot.title = element_text(size=size_title),
plot.caption = element_text(hjust = 0, face= "italic",size=9))+
ggtitle("b) Respiratory Consultations")+
scale_x_date(date_breaks = "3 months",date_labels="%b %y",limits=c(as.Date("2015-01-05"),as.Date("2019-12-23")))
rate_trauma_plot2_red<-
plot(impact3d_ratio, "pointwise")$data %>%
dplyr::left_join(data15a64_rn[,c("date","rn")], by=c("time"="rn")) %>%
ggplot(aes(x=date))+
geom_line(aes(y=mean), size=line_size, linetype="longdash", color="darkblue") +
geom_line(aes(y=response), size=line_size, color="black")+
geom_ribbon(aes(ymin=lower,ymax=upper),fill="steelblue", alpha = 0.35)+
theme_sjplot()+
geom_hline(yintercept = 0, col = "black", lty = 3, size=1)+
geom_vline(xintercept = length(data15a64_rn$did[data15a64_rn$did==0])+1, col = "gray50", lty = 2, size=1)+
theme(axis.text.x = element_text(angle = angle_x, vjust = 0.5, hjust=.5,size=height_x),
# theme(axis.text.x = element_blank(),#element_text(angle = 90, vjust = 0.5, hjust=.5,size=9),
axis.text.y = element_text(size=height_y),
axis.title.y = element_blank(),
axis.title.x = element_blank(),
plot.title = element_text(size=size_title),
plot.caption = element_text(hjust = 0, face= "italic",size=9))+
ggtitle("e) Trauma Hospitalizations per 1,000 consultations")+
scale_x_date(date_breaks = "3 months",date_labels="%b %y",limits=c(as.Date("2015-01-05"),as.Date("2019-12-23")))
rate_resp_plot2_red<-
plot(impact3d_ratio_resp, "pointwise")$data %>%
dplyr::left_join(data15a64_rn[,c("date","rn")], by=c("time"="rn")) %>%
ggplot(aes(x=date))+
geom_line(aes(y=mean), size=line_size, linetype="longdash", color="darkblue") +
geom_line(aes(y=response), size=line_size, color="black")+
geom_ribbon(aes(ymin=lower,ymax=upper),fill="steelblue", alpha = 0.35)+
theme_sjplot()+
geom_hline(yintercept = 0, col = "black", lty = 3, size=1)+
geom_vline(xintercept = length(data15a64_rn$did[data15a64_rn$did==0])+1, col = "gray50", lty = 2, size=1)+
theme(axis.text.x = element_text(angle = angle_x, vjust = 0.5, hjust=.5,size=height_x),
# theme(axis.text.x = element_blank(),#element_text(angle = 90, vjust = 0.5, hjust=.5,size=9),
axis.text.y = element_text(size=height_y),
axis.title.y = element_blank(),
axis.title.x = element_blank(),
plot.title = element_text(size=size_title),
plot.caption = element_text(hjust = 0, face= "italic",size=9))+
ggtitle("f) Respiratory Hospitalizations per 1,000 consultations")+
scale_x_date(date_breaks = "3 months",date_labels="%b %y",limits=c(as.Date("2015-01-05"),as.Date("2019-12-23")))
vector_breaks_plot2b_red<- seq(238,nrow(data15a64_rn),3)
Sys.setlocale(category = "LC_ALL", locale = "english")## [1] "LC_COLLATE=English_United States.1252;LC_CTYPE=English_United States.1252;LC_MONETARY=English_United States.1252;LC_NUMERIC=C;LC_TIME=English_United States.1252"vector_dates_plot2b_red<-as.character(format(as.Date(unlist(data15a64_rn[c(seq(238,nrow(data15a64_rn),3)),"date"])),"%b %d"))
plot2b_red<-
plot_grid(
cons_trauma_plot2_red+scale_y_continuous(breaks=seq(-800,250,200), limits = c(-800,250)),
cons_resp_plot2_red+scale_y_continuous(breaks=seq(-800,250,200), limits = c(-800,250)),
hosp_trauma_plot2_red+scale_y_continuous(breaks=seq(-75,75,25), limits = c(-75,75)),
hosp_resp_plot2_red+scale_y_continuous(breaks=seq(-75,75,25), limits = c(-75,75)),
rate_trauma_plot2_red+scale_y_continuous(breaks=seq(-200,400,200), limits = c(-200,400)),
rate_resp_plot2_red+scale_y_continuous(breaks=seq(-200,400,200), limits = c(-200,400)),
#get_legend(t14_trend_leg + theme(legend.position="right",
# legend.text = element_text(size = 15))),
nrow = 3, rel_widths = c(.95, .95), rel_heights = c(rep(1,2),1.3), scale = c(.95, .95, .95,.95,.95,.95)
)+
draw_label("Weekly difference change in consultations/hospitalizations",x=0.005, y=.5, angle=90, size = 17)
#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:
#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:
#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:
ggdraw(plot2b_red)+
theme(plot.background = element_rect(fill=NA, color = NA))
Figure 2. Estimated Trends of the Outcomes
jpeg("_Fig1_final.jpg", height = 14, width = 23, units = 'in', res = 600)
ggdraw(plot2b_red)+
theme(plot.background = element_rect(fill=NA, color = NA))
dev.off()## png
## 2cairo_ps("_Fig1_final.eps",
height = 14, width = 23,
fallback_resolution = 600,
family= "Times")
#horizontal=T,
#paper= "letter",
#pagecentre=F)
ggdraw(plot2b_red)+
theme(plot.background = element_rect(fill=NA, color = NA))
dev.off()## png
## 2 ggsave("_Fig1_final.ps", plot =
ggdraw(plot2b_red)+
theme(plot.background = element_rect(fill=NA, color = NA)),
path = NULL,
scale = 1, width = 23, height = 14, units = "in",
device=cairo_ps,
fallback_resolution = 600,
limitsize = F)#horizontal
height_y<-16
height_x<-16
size_title<-18
line_size<-.8
angle_x <- 45
hosp_trauma_plot3_red<-
plot(impact3d_hosp_trauma, "cumulative")$data %>%
dplyr::left_join(data15a64_rn[,c("date","rn")], by=c("time"="rn")) %>%
ggplot(aes(x=date))+
geom_line(aes(y=mean), size=line_size, linetype="longdash", color="darkblue") +
geom_line(aes(y=response), size=line_size, color="black")+
geom_ribbon(aes(ymin=lower,ymax=upper),fill="steelblue", alpha = 0.35)+
theme_sjplot()+
geom_hline(yintercept = 0, col = "black", lty = 3, size=1)+
geom_vline(xintercept = length(data15a64_rn$did[data15a64_rn$did==0])+1, col = "gray50", lty = 2, size=1)+
theme(axis.text.x = element_blank(),
# theme(axis.text.x = element_blank(),#element_text(angle = 90, vjust = 0.5, hjust=.5,size=9),
axis.text.y = element_text(size=height_y),
axis.title.y = element_blank(),
axis.title.x = element_blank(),
plot.title = element_text(size=size_title),
plot.caption = element_text(hjust = 0, face= "italic",size=9))+
ggtitle("c) Trauma Hospitalizations")+
scale_x_date(date_breaks = "3 months",date_labels="%b %y",limits=c(as.Date("2015-01-05"),as.Date("2019-12-23")))
hosp_resp_plot3_red<-
plot(impact3d1_hosp_resp, "cumulative")$data %>%
dplyr::left_join(data15a64_rn[,c("date","rn")], by=c("time"="rn")) %>%
ggplot(aes(x=date))+
geom_line(aes(y=mean), size=line_size, linetype="longdash", color="darkblue") +
geom_line(aes(y=response), size=line_size, color="black")+
geom_ribbon(aes(ymin=lower,ymax=upper),fill="steelblue", alpha = 0.35)+
theme_sjplot()+
geom_hline(yintercept = 0, col = "black", lty = 3, size=1)+
geom_vline(xintercept = length(data15a64_rn$did[data15a64_rn$did==0])+1, col = "gray50", lty = 2, size=1)+
theme(axis.text.x = element_blank(),
# theme(axis.text.x = element_blank(),#element_text(angle = 90, vjust = 0.5, hjust=.5,size=9),
axis.text.y = element_text(size=height_y),
axis.title.y = element_blank(),
axis.title.x = element_blank(),
plot.title = element_text(size=size_title),
plot.caption = element_text(hjust = 0, face= "italic",size=9))+
ggtitle("d) Respiratory Hospitalizations")+
scale_x_date(date_breaks = "3 months",date_labels="%b %y",limits=c(as.Date("2015-01-05"),as.Date("2019-12-23")))
cons_trauma_plot3_red<-
plot(impact3d1_cons_trauma, "cumulative")$data %>%
dplyr::left_join(data15a64_rn[,c("date","rn")], by=c("time"="rn")) %>%
ggplot(aes(x=date))+
geom_line(aes(y=mean), size=line_size, linetype="longdash", color="darkblue") +
geom_line(aes(y=response), size=line_size, color="black")+
geom_ribbon(aes(ymin=lower,ymax=upper),fill="steelblue", alpha = 0.35)+
theme_sjplot()+
geom_hline(yintercept = 0, col = "black", lty = 3, size=1)+
geom_vline(xintercept = length(data15a64_rn$did[data15a64_rn$did==0])+1, col = "gray50", lty = 2, size=1)+
theme(axis.text.x = element_blank(),
# theme(axis.text.x = element_blank(),#element_text(angle = 90, vjust = 0.5, hjust=.5,size=9),
axis.text.y = element_text(size=height_y),
axis.title.y = element_blank(),
axis.title.x = element_blank(),
plot.title = element_text(size=size_title),
plot.caption = element_text(hjust = 0, face= "italic",size=9))+
ggtitle("a) Trauma Consultations")+
scale_x_date(date_breaks = "3 months",date_labels="%b %y",limits=c(as.Date("2015-01-05"),as.Date("2019-12-23")))
cons_resp_plot3_red<-
plot(impact3d_cons_resp, "cumulative")$data %>%
dplyr::left_join(data15a64_rn[,c("date","rn")], by=c("time"="rn")) %>%
ggplot(aes(x=date))+
geom_line(aes(y=mean), size=line_size, linetype="longdash", color="darkblue") +
geom_line(aes(y=response), size=line_size, color="black")+
geom_ribbon(aes(ymin=lower,ymax=upper),fill="steelblue", alpha = 0.35)+
theme_sjplot()+
geom_hline(yintercept = 0, col = "black", lty = 3, size=1)+
geom_vline(xintercept = length(data15a64_rn$did[data15a64_rn$did==0])+1, col = "gray50", lty = 2, size=1)+
theme(axis.text.x = element_blank(),
# theme(axis.text.x = element_blank(),#element_text(angle = 90, vjust = 0.5, hjust=.5,size=9),
axis.text.y = element_text(size=height_y),
axis.title.y = element_blank(),
axis.title.x = element_blank(),
plot.title = element_text(size=size_title),
plot.caption = element_text(hjust = 0, face= "italic",size=9))+
ggtitle("b) Respiratory Consultations")+
scale_x_date(date_breaks = "3 months",date_labels="%b %y",limits=c(as.Date("2015-01-05"),as.Date("2019-12-23")))
rate_trauma_plot3_red<-
plot(impact3d_ratio, "cumulative")$data %>%
dplyr::left_join(data15a64_rn[,c("date","rn")], by=c("time"="rn")) %>%
ggplot(aes(x=date))+
geom_line(aes(y=mean), size=line_size, linetype="longdash", color="darkblue") +
geom_line(aes(y=response), size=line_size, color="black")+
geom_ribbon(aes(ymin=lower,ymax=upper),fill="steelblue", alpha = 0.35)+
theme_sjplot()+
geom_hline(yintercept = 0, col = "black", lty = 3, size=1)+
geom_vline(xintercept = length(data15a64_rn$did[data15a64_rn$did==0])+1, col = "gray50", lty = 2, size=1)+
theme(axis.text.x = element_text(angle = angle_x, vjust = 0.5, hjust=.5,size=height_x),
# theme(axis.text.x = element_blank(),#element_text(angle = 90, vjust = 0.5, hjust=.5,size=9),
axis.text.y = element_text(size=height_y),
axis.title.y = element_blank(),
axis.title.x = element_blank(),
plot.title = element_text(size=size_title),
plot.caption = element_text(hjust = 0, face= "italic",size=9))+
ggtitle("e) Trauma Hospitalizations per 1,000 consultations")+
scale_x_date(date_breaks = "3 months",date_labels="%b %y",limits=c(as.Date("2015-01-05"),as.Date("2019-12-23")))
rate_resp_plot3_red<-
plot(impact3d_ratio_resp, "cumulative")$data %>%
dplyr::left_join(data15a64_rn[,c("date","rn")], by=c("time"="rn")) %>%
ggplot(aes(x=date))+
geom_line(aes(y=mean), size=line_size, linetype="longdash", color="darkblue") +
geom_line(aes(y=response), size=line_size, color="black")+
geom_ribbon(aes(ymin=lower,ymax=upper),fill="steelblue", alpha = 0.35)+
theme_sjplot()+
geom_hline(yintercept = 0, col = "black", lty = 3, size=1)+
geom_vline(xintercept = length(data15a64_rn$did[data15a64_rn$did==0])+1, col = "gray50", lty = 2, size=1)+
theme(axis.text.x = element_text(angle = angle_x, vjust = 0.5, hjust=.5,size=height_x),
# theme(axis.text.x = element_blank(),#element_text(angle = 90, vjust = 0.5, hjust=.5,size=9),
axis.text.y = element_text(size=height_y),
axis.title.y = element_blank(),
axis.title.x = element_blank(),
plot.title = element_text(size=size_title),
plot.caption = element_text(hjust = 0, face= "italic",size=9))+
ggtitle("f) Respiratory Hospitalizations per 1,000 consultations")+
scale_x_date(date_breaks = "3 months",date_labels="%b %y",limits=c(as.Date("2015-01-05"),as.Date("2019-12-23")))
vector_breaks_plot3b_red<- seq(238,nrow(data15a64_rn),3)
Sys.setlocale(category = "LC_ALL", locale = "english")## [1] "LC_COLLATE=English_United States.1252;LC_CTYPE=English_United States.1252;LC_MONETARY=English_United States.1252;LC_NUMERIC=C;LC_TIME=English_United States.1252"vector_dates_plot3b_red<-as.character(format(as.Date(unlist(data15a64_rn[c(seq(238,nrow(data15a64_rn),3)),"date"])),"%b %d"))
plot3b_red<-
plot_grid(
cons_trauma_plot3_red+scale_y_continuous(breaks=seq(-4000,1500,1000), limits = c(-4000,1500)),
cons_resp_plot3_red+scale_y_continuous(breaks=seq(-4000,1500,1000), limits = c(-4000,1500)),
hosp_trauma_plot3_red+scale_y_continuous(breaks=seq(-200,200,100), limits = c(-200,200)),
hosp_resp_plot3_red+scale_y_continuous(breaks=seq(-200,200,100), limits = c(-200,200)),
rate_trauma_plot3_red+scale_y_continuous(breaks=seq(-100,1250,250), limits = c(-100,1250)),
rate_resp_plot3_red+scale_y_continuous(breaks=seq(-100,1250,250), limits = c(-100,1250)),
#get_legend(t14_trend_leg + theme(legend.position="right",
# legend.text = element_text(size = 15))),
nrow = 3, rel_widths = c(1, 1), rel_heights = c(rep(1,2),1.25), scale = c(.95, .95, .95,.95,.95,.95))+
draw_label("Weekly difference change in consultations/hospitalizations",x=0.005, y=.5, angle=90, size = 17)
#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:
#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:
#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:
ggdraw(plot3b_red)+
theme(plot.background = element_rect(fill=NA, color = NA))
Figure 3. Estimated Trends of the Outcomes
jpeg("_Fig2_final.jpg", height = 14, width = 23, units = 'in', res = 600)
ggdraw(plot3b_red)+
theme(plot.background = element_rect(fill=NA, color = NA))
dev.off()## png
## 2cairo_ps("_Fig2_final.eps",
height = 14, width = 23,
fallback_resolution = 600,
family= "Times")
#horizontal=T,
#paper= "letter",
#pagecentre=F)
ggdraw(plot3b_red)+
theme(plot.background = element_rect(fill=NA, color = NA))
dev.off()## png
## 2 ggsave("_Fig2_final.ps", plot =
ggdraw(plot3b_red)+
theme(plot.background = element_rect(fill=NA, color = NA)),
path = NULL,
scale = 1, width = 23, height = 14, units = "in",
device=cairo_ps,
fallback_resolution = 600,
limitsize = F)0.2 Consolidated Plot
# dejar una sola figura que contenga 6 lÃneas (una para cada outcome) con el cambio porcentual (para que todos los outcomes estén en la misma unidad de medida).
#_#_#_#_#_# CONS TRAUMA #_#_#_#_#_#
releff_cons_trauma<-
cbind(plot(impact3d1_cons_trauma, "original")$data,
data15a64_rn[,c("date","did")])%>%
dplyr::mutate_at(.vars=vars(response,mean,lower,upper),
.funs = funs(`zw`=scale(.))) %>%
dplyr::mutate(diff=response-mean) %>%
dplyr::mutate(rel_diff=diff/mean) %>%
dplyr::mutate(diff_zw=response_zw-mean_zw) %>%
#ggplot(aes(y=rel_diff, x=time))+
# geom_line()+
# geom_vline(xintercept = 253, color="red")+
# theme_bw()
dplyr::select(date,time,rel_diff,diff_zw) #%>%
#tail(20)
#ggplot(aes(y=rel_diff, x=date))+
#geom_line()
# impact3d1_cons_trauma
# Relative effect (s.d.) -14% (13%) -14% (13%)
# 95% CI [-40%, 12%] [-40%, 12%]
#_#_#_#_#_# CONS RESP #_#_#_#_#_#
releff_cons_resp<-
cbind(plot(impact3d_cons_resp, "original")$data,
data15a64_rn[,c("date","did")])%>%
dplyr::mutate_at(.vars=vars(response,mean,lower,upper),
.funs = funs(`zw`=scale(.))) %>%
dplyr::mutate(diff=response-mean) %>%
dplyr::mutate(rel_diff=diff/mean) %>%
dplyr::mutate(diff_zw=response_zw-mean_zw) %>%
dplyr::select(date,time,rel_diff,diff_zw) #%>%
#tail(20)
#ggplot(aes(y=rel_diff, x=date))+
#geom_line()
# impact3d_cons_resp
# Relative effect (s.d.) -30% (30%) -30% (30%)
# 95% CI [-89%, 30%] [-89%, 30%]
#_#_#_#_#_# HOSP TRAUMA #_#_#_#_#_#
releff_hosp_trauma<-
cbind(plot(impact3d_hosp_trauma, "original")$data,
data15a64_rn[,c("date","did")])%>%
dplyr::mutate_at(.vars=vars(response,mean,lower,upper),
.funs = funs(`zw`=scale(.))) %>%
dplyr::mutate(diff=response-mean) %>%
dplyr::mutate(rel_diff=diff/mean) %>%
dplyr::mutate(diff_zw=response_zw-mean_zw) %>%
dplyr::select(date,time,rel_diff,diff_zw) #%>%
#tail(20)
#ggplot(aes(y=rel_diff, x=date))+
#geom_line()
# impact3d_cons_resp
# Relative effect (s.d.) 15% (5.6%) 15% (5.6%)
# 95% CI [4%, 26%] [4%, 26%]
#_#_#_#_#_# HOSP RESP #_#_#_#_#_#
releff_hosp_resp<-
cbind(plot(impact3d1_hosp_resp, "original")$data,
data15a64_rn[,c("date","did")])%>%
dplyr::mutate_at(.vars=vars(response,mean,lower,upper),
.funs = funs(`zw`=scale(.))) %>%
dplyr::mutate(diff=response-mean) %>%
dplyr::mutate(rel_diff=diff/mean) %>%
dplyr::mutate(diff_zw=response_zw-mean_zw) %>%
dplyr::select(date,time,rel_diff,diff_zw) #%>%
#tail(20)
#ggplot(aes(y=rel_diff, x=date))+
#geom_line()
# impact3d_cons_resp
# Relative effect (s.d.) -6.8% (19%) -6.8% (19%)
# 95% CI [-44%, 31%] [-44%, 31%]
#_#_#_#_#_# RATIO TRAUMA #_#_#_#_#_#
releff_ratio_trauma<-
cbind(plot(impact3d_ratio, "original")$data,
data15a64_rn[,c("date","did")])%>%
dplyr::mutate_at(.vars=vars(response,mean,lower,upper),
.funs = funs(`zw`=scale(.))) %>%
dplyr::mutate(diff=response-mean) %>%
dplyr::mutate(rel_diff=diff/mean) %>%
dplyr::mutate(diff_zw=response_zw-mean_zw) %>%
dplyr::select(date,time,rel_diff,diff_zw) #%>%
#tail(20)
#ggplot(aes(y=rel_diff, x=date))+
#geom_line()
# impact3d_cons_resp
# Relative effect (s.d.) 40% (14%) 40% (14%)
# 95% CI [13%, 68%] [13%, 68%]
#_#_#_#_#_# RATIO RESP #_#_#_#_#_#
releff_ratio_resp<-
cbind(plot(impact3d_ratio_resp, "original")$data,
data15a64_rn[,c("date","did")])%>%
dplyr::mutate_at(.vars=vars(response,mean,lower,upper),
.funs = funs(`zw`=scale(.))) %>%
dplyr::mutate(diff=response-mean) %>%
dplyr::mutate(rel_diff=diff/mean) %>%
dplyr::mutate(diff_zw=response_zw-mean_zw) %>%
dplyr::select(date,time,rel_diff,diff_zw) #%>%
#tail(20)
#ggplot(aes(y=rel_diff, x=date))+
#geom_line()
# impact3d_cons_resp
# Relative effect (s.d.) 59% (15%) 59% (15%)
# 95% CI [29%, 88%] [29%, 88%]#horizontal
height_y<-18
height_x<-18
size_title<-18
line_size<-.8
angle_x <- 45
Sys.setlocale(category = "LC_ALL", locale = "english")## [1] "LC_COLLATE=English_United States.1252;LC_CTYPE=English_United States.1252;LC_MONETARY=English_United States.1252;LC_NUMERIC=C;LC_TIME=English_United States.1252"manualcolors<-c('indianred1','cornflowerblue', 'gray50', 'darkolivegreen4', 'slateblue2', 'firebrick4', 'goldenrod4')
releff<-
releff_cons_trauma[,c("date","rel_diff")] %>%
dplyr::left_join(releff_cons_resp[,c("date","rel_diff")],by="date") %>%
dplyr::left_join(releff_hosp_trauma[,c("date","rel_diff")],by="date") %>%
dplyr::left_join(releff_hosp_resp[,c("date","rel_diff")],by="date") %>%
dplyr::left_join(releff_ratio_trauma[,c("date","rel_diff")],by="date") %>%
dplyr::left_join(releff_ratio_resp[,c("date","rel_diff")],by="date")
colnames(releff)<-c("date","Cons_Trauma","Cons_Resp","Hosp_Trauma","Hosp_Resp","Ratio_Trauma","Ratio_Resp")
releff_melt<-
releff %>%
melt(id.vars=c("date")) %>%
dplyr::mutate(variable=
dplyr::case_when(variable=="Cons_Trauma"~"Trauma Consultations",
variable=="Cons_Resp"~"Respiratory Consultations",
variable=="Hosp_Trauma"~"Trauma Hospitalizations",
variable=="Hosp_Resp"~"Respiratory Hospitalizations",
variable=="Ratio_Trauma"~"Trauma Hospitalizations per 1,000 consultations",
variable=="Ratio_Resp"~"Respiratory Hospitalizations per 1,000 consultations")) %>%
dplyr::filter(date>="2019-08-05")
releff_plot<-
releff_melt%>%
dplyr::mutate(variable=factor(variable)) %>%
ggplot(aes(x=date, y=value, color=variable))+ #, shape=variable
geom_line(size=line_size)+
#scale_shape_manual(name= "Outcome", values=c(1:6))+
#scale_linetype_manual(name= "",
# values= c("twodash","solid","dotted","dashed","twodash","solid"))+
scale_color_manual(name= "", values=manualcolors)+
theme_sjplot()+
theme(legend.position = "bottom",
legend.text= element_text(size=18),
axis.text.y= element_text(size=height_y),
axis.title.y= element_text(size=height_y),
axis.text.x = element_text(angle = angle_x, vjust = 0.5, hjust=.5,size=height_x))+
labs(y="Percent change", x="")+
scale_x_date(breaks = scales::date_breaks("4 weeks"), labels = scales::date_format("%d-%b"))+
scale_y_continuous(labels = scales::percent_format(accuracy = 1))+
geom_vline(xintercept = as.Date("2019-10-21"), col = "darkred", lty = 2, size=1)+
geom_hline(yintercept = 0, col = "black", lty = 3, size=1)
releff_plot
Figure 4. Relative Effects
cairo_ps("_Fig1_final_AJPH.eps",
height = 14, width = 23,
fallback_resolution = 600,
family= "Times")
#horizontal=T,
#paper= "letter",
#pagecentre=F)
releff_plot+
theme(plot.background = element_rect(fill=NA, color = NA))
dev.off()## png
## 2 ggsave("_Fig1_final_AJPH.ps", plot =
releff_plot+
theme(plot.background = element_rect(fill=NA, color = NA)),
path = NULL,
scale = 1, width = 23, height = 14, units = "in",
device=cairo_ps,
fallback_resolution = 600,
limitsize = F)
jpeg("_Fig1_AJPH.jpg", height = 14, width = 23, units = 'in', res = 600)
releff_plot+
theme(plot.background = element_rect(fill=NA, color = NA))
dev.off()## png
## 2